# Code to derive first differences from MCMCglmm model


library(postMCMCglmm)


###################################################################################################################
# PREDICTED PROBABILITIES [3 CAT]
###################################################################################################################

predict2.MCMCglmm.3cat <- predict2

ests.MCMCglmm.3cat <- function(res){
	ans1 <- c(mean(res[[1]]), quantile(res[[1]], probs=c(.025, .975)))
	ans2 <- c(mean(res[[2]]), quantile(res[[2]], probs=c(.025, .975)))
	ans3 <- c(mean(res[[3]]), quantile(res[[3]], probs=c(.025, .975)))

	ans <- as.data.frame(rbind(ans1, ans2, ans3))
	names(ans) <- c("mean", "lwr", "upr")
	row.names(ans) <- c("-1", "0", "+1")
	ans	
}


###################################################################################################################
# FIRST DIFFERENCES 
###################################################################################################################

fd.MCMCglmm.3cat <- function(object, X, X1, Z, use = c("all", "mean"), type = c("lp", "response"), ...) {

# SET-UP
  use <- match.arg(use)
  type <- match.arg(type)

  if (!is.null(object$X) && missing(X)) {
    X <- object$X
  }
  if (!is.null(object$Z) && missing(Z)) {
    Z <- object$Z
  }

  if (missing(X)) X <- NULL
  if (missing(Z)) Z <- NULL
  Xb <- Zu <- 0L

  b <- fixef(object, use = use)
  u <- ranef(object, use = use)
  stddev <- sqrt(unique(object$error.term) + 1)



# PRED 1

  if (!is.null(X)) Xb <- X %*% b
  if (!is.null(Z)) Zu <-  Z %*% u
  res <- t(as.matrix(Xb + Zu))
  
  dimnames(res) <- list(switch(use,
    all = paste0("Rep.", 1:nrow(res)), mean = NULL), NULL)

      # cut points
      CP <- object$CP
      if (use == "mean") CP <- colMeans(CP)
      CP <- as.list(as.data.frame(CP))
      CP <- c(0, lapply(CP, function(x) {
        do.call("cbind", rep(list(x), dim(res)[2L]))
      }))
      # difference between cuts and predicted plus lower and upper bounds
      for (i in seq_along(CP)) {
        CP[[i]] <- pnorm(CP[[i]] - res, 0, stddev)
      }
      CP <- c(0, CP, 1)
      q <- vector("list", length(CP) - 2)
      for (i in 2:(length(CP) - 1)) {
        q[[i - 1]] <- mcmc(CP[[i + 1]] - CP[[i]])
      }
      q1 <- c(list(mcmc(1 - Reduce(`+`, q[1:(i - 1)]))), q)
      q1 <- list(q1[[1]], q1[[2]], q1[[3]])
      
      
# PRED 2

  if (!is.null(X1)) Xb <- X1 %*% b
  if (!is.null(Z)) Zu <-  Z %*% u
  res <- t(as.matrix(Xb + Zu))
  
  dimnames(res) <- list(switch(use,
    all = paste0("Rep.", 1:nrow(res)), mean = NULL), NULL)

      # cut points
      CP <- object$CP
      if (use == "mean") CP <- colMeans(CP)
      CP <- as.list(as.data.frame(CP))
      CP <- c(0, lapply(CP, function(x) {
        do.call("cbind", rep(list(x), dim(res)[2L]))
      }))
      # difference between cuts and predicted plus lower and upper bounds
      for (i in seq_along(CP)) {
        CP[[i]] <- pnorm(CP[[i]] - res, 0, stddev)
      }
      CP <- c(0, CP, 1)
      q <- vector("list", length(CP) - 2)
      for (i in 2:(length(CP) - 1)) {
        q[[i - 1]] <- mcmc(CP[[i + 1]] - CP[[i]])
      }
      q2 <- c(list(mcmc(1 - Reduce(`+`, q[1:(i - 1)]))), q)
      q2 <- list(q2[[1]], q2[[2]], q2[[3]])      
      
      
      
# FDs

	fd1 <- q1[[1]]-q2[[1]]
	fd2 <- q1[[2]]-q2[[2]]
	fd3 <- q1[[3]]-q2[[3]]
	res <- list(fd1,fd2,fd3)
    res <- as.mcmc(res)
    class(res) <- c("mcmc", "MCMCglmmPredictedLP")
 	return(res)
}
   
